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SUMMARY 

This paper presents the results of a numerical simulation of the time- 
averaged Invlscld flow field through the blade rows of a multiblade row 
turboprop configuration. The governing equations are outlined along with a 
discussion of the solution procedure and coding strategy. Numerical results 
obtained from a simulation of the flow field through a modern high-speed 
turboprop will be shown. 


INTRODUCTION 

Solving the flow field of a multiblade row machine Is very difficult con- 
sidering the time and length scales Involved. The numerical simulation of a 
general configuration remains formidable even for a machine much larger and 
faster than today's Class VI computers. However, by mathematically modeling 
the flow field In the spirit of Reynolds-averaged modeling of turbulent flows, 
much of the physics relevant to design can be deduced (ref. 1). The objective 
of this paper Is to describe a procedure for simulating the Invlscld, time- 
average flow through a multiblade row geometry and present results for a 
counter-rotating propeller configuration designed to operate at transonic 
speeds. 

Three-dimensional, Invlscld codes have been developed for Isolated propel- 
lers using a variety of algorithms. Bober, et al. (ref. 2) and Barton, et al. 
(ref. 3) used the Beam and Warming algorithm and obtained good comparisons to 
experimental results. Clark (ref. 4) used Denton's finite volume code modified 
for propellers to obtain solutions for acoustical analyses. Holmes and Tong 
(ref. 5) applied Jameson's Runge-Kutta procedure (ref. 6) formulated In terms 
of cartesian velocity components to a turbine, compressor, and propeller blade 
row. Celestina and Adamczyk (ref. 7) presented results applying Jameson's 
technique formulated In terms of cylindrical velocity components to a turbine 
and propeller blade row. 

A procedure for extending Isolated blade row analyses to multiblade row 
configurations was suggested by Denton (ref. 8). His method as applied to a 
stage Involved circumferentially averaging the flow properties at a given axial 
location between the two blade rows. By doing so, the downstream boundary 



condition to the first blade row and upstream boundary condition to the second 
blade row are circumferentially uniform. 

Adamczyk (ref. 1) developed a first principles procedure for analyzing 
multiblade row flows In which a sequence of averaging operators are used to 
derive a set of equations that describes an "averaged" three-dimensional rep- 
resentation of the flow field through each blade row of a multiblade row 
machine. This flow field Is steady In time In a reference frame fixed to each 
blade row and spatially periodic from passage to passage. The field equations 
associated with this model are referred to as the average-passage equation 
system. The present work Is therefore unique In that It outlines a procedure 
to solve a set of equations which govern a three-dimensional multiblade row 
flow field that Is directly traceable to the Navler-Stokes equations.. However, 
this average-passage description has associated with It a well known diffi- 
culty. Averaged equations of motion always lead to situations In which there 
are more unknowns than equations - this Is the "closure" problem which requires 
assumptions or empirical Information to make the number of equations equal to 
the number of unknowns. The closure problem Is addressed for the Invlscld form 
of the equations by Adamczyk In reference 9. 

The procedure for solving the average-passage equation system requires a 
mesh for each blade row. The axial and radial location of each grid point Is 
Identical for each mesh, but varies In the circumferential direction due to 
relative blade row locations and unequal blade numbers. A solution Is gener- 
ated for each blade row on Its own mesh with the effect of the neighboring 
blade row contained In source terms In the governing equations. The source 
terms are sequentially updated using Information provided by the neighboring 
blade row simulation. 


GOVERNING EQUATIONS 

The three-dimensional average-passage equation system for simulating the 
flow through multiblade rows can be written In cylindrical (r,e,z) coordinates 
as 

(Mi) t + L(xu.) + J*XS dVol = J XK dVol (1) 

The vector u contains the flow variables density, axial and radial momenta, 
angular momentum and total Internal energy and X Is the neighboring blade row 
blockage factor. The value of this parameter ranges between zero and unity, 
unity being the value associated with zero blade thickness. (See ref. 1 for 
details.) The operator l(Xu) balances the mass, axial and radial momenta, 
angular momentum, and energy through a control volume, JXK dVol Is a source 
term due to the cylindrical coordinate system, and JS dVol contains the body 
forces, energy sources, momentum, and energy temporal mixing correlations 
associated with the neighboring blade row(s). The details for computing these 
terms are given In a companion paper (ref. 9). The vector u and the oper- 
ator L( Xu) are defined by the following expressions: 
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In the above equations p represents the density, v r , v 0 , and v z are the 
radial, tangential, and axial absolute velocities and p Is the pressure. 
From the equation of state the total Internal energy Is related to the 
pressure as follows 
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and the total enthalpy Is related to e 0 , p and by the equation 





In the above equations all lengths are nondlmenslonallzed by the diameter of 
the largest blade row. The velocity components are nondlmenslonallzed by the 
reference speed of sound, a re f/v/y, pressure and density by their respective 
reference values, and y Is the ratio of specific heats. 

For rotating flows, the absolute (fixed) reference frame Is transformed 
to the relative (rotating) frame by the transformation 


e . . . = e , + ftt (3) 

absolute relative 

where Q is the rotational wheel speed (positive with e) . Introducing 
equation (3) Into equation (1) yields 

L(Xu) XF • dA p + \(G - rnu) • dA Q + \H • dA z (4) 

The term (G - r«u) • d^ represents the relative flux of u In the tangential 
direction. 


Equation (1) Is discretized In space for a cell volume (fig. 1) by 
approximating the surface Integrals by the mid point rule. The result Is a 
system of ordinary differential equations of the form 
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(Xu) ^ |xF dA r + X( G - rQu) dA^ + \HdA z J + J \S dVol = j \Kd Vol (5) 

The surface areas, dA r ,dAe,dA z are calculated using the cross product of the 
diagonals of a cell face and the volume Is determined using the formula des- 
cribed by Holmes and Tong (ref. 5). Since all the flow quantities are cell- 
centered, a simple averaging procedure Is used to determine the value of a 
variable at any surface, excluding solid boundaries. This Is equivalent to 
second-order accurate central differencing for a uniform mesh. 


Runge-Kutta Integration 


To advance the equations In time, a four-stage Runge-Kutta scheme Is used. 
The scheme employed has been patterned after the work of Jameson, et al. 

(ref. 6). Given Information at time level n, the steps to advance to the 
next level , n + 1 , are 
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where a = 1/4, 3 = 1/3, y = 1/2 and D(u) Is the dissipation operator. The 
maximum permissible time step for this scheme Is restricted by the CFL stabil- 
ity limit. Jameson has determined the limit for the above four-stage scheme 
to be 2/7 based on a one-dimensional model problem. To enhance the conver- 
gence rate of this scheme, a local time step Is chosen based on the maximum CFL 
number commensurate with stability. An advantage of using the present Runge- 
Kutta scheme Is that It minimizes storage requirements. 


Artificial Dissipation 


To suppress odd-even point decoupling In the solution, dissipative terms 
are added to the equations. Jameson (ref. 6), via numerical experiments, 
developed a blend of second and fourth difference smoothing operators. The 
operator D(u) In equation (6) can be decomposed Into three spatial operators 

D(u) = ( D + D + D ) (u) 
r v7 z 


such that the dissipation In each direction can be evaluated separately. The 
dissipation In the axial direction, D z (u), Is expressed as follows 


D z (u) 
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and the coefficients c and e are evaluated as follows: 

c 1+l/2,j ,k = K2max(v 1+l,J,k, w 1,j,k ) , 

°’ ,c Hl/2,j,k " e 1+l/2,J , 

k 2 and k 4 are constants typically set at 1/2 and 1/64, respectively. 

To capture shocks sharply and retain second-order accuracy away from shocks, 
Jameson defined the coefficient wi,j,ic as 
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The variable j k Is proportional to the square of the mesh spacing In 
smooth regions of flow and linear In mesh spacing In regions of large pressure 
gradients. The pressure sensitive switch was applied only In the axial and 
circumferential directions. For the radial direction, It was found that set- 
ting "1,1, k to a constant (0.05) and e 4 to zero enhanced the stability 
of the scneme. ' 

■ A 

Boundary Conditions 

The mathematical form of the field equations solved In the present work 
closely resembles that of the Euler equations. Therefore, the known mathemati- 
cal properties of the Euler equations will be used as a guide to develop the 
boundary conditions of the present system of equations. It Is assumed that the 
absolute flow field approaching and leaving the propeller Is subsonic. This 
Implies that four conditions must be specified at the upstream boundary and 
only one at the downstream boundary. 


The axial velocity component and the flow properties at the upstream 
boundary are updated based on a local unsteady one-dimensional flow model In 
which the entropy Is assumed constant with time and uniform In space. The 
equations associated with this model are 
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In which C + and C~ are the well-known Rlemann Invariants. They are related 
to the axial velocity, v z , and speed of sound, a, by the equations: 



The Invariant C + Is associated with Information coming from outside the 
computational domain and thus It Is specified based on the far field flow 
conditions. The C - Invariant Is updated by solving equation (10b) In time 
using the Runge-Kutta Integration procedure outlined earlier. The axial 
derivative Is approximated by a backward difference operator. The Rlemann 
Invariants determine the speed of sound and the axial velocity component. The 
pressure, density, and temperature are updated based on the known value of the 
Incoming entropy. The value of the velocity components parallel to the 
upstream boundary are assumed known. 

At the downstream boundary, simple radial equilibrium Is enforced. 
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The pressure Is specified at the free-stream farfleld boundary and equa- 
tion (11) Is integrated radially using the trapezoidal rule toward the 
splnner-nacel le. The remaining flow variables are extrapolated from the 
Interior. 


Boundary conditions at the farfleld boundary are derived based on a one- 
dimensional unsteady flow model Identical to the one employed at the upstream 
boundary. In this model, the axial velocity component Is replaced by the 
radial component. The value of C+ Is fixed by the farfleld condition and 
C~ Is extrapolated from the Interior. 

At periodic flow boundaries we require the flow to exhibit a spatial 
periodicity equal to the pitch of the blade row. Thus, any Information 
required from a cell which lies adjacent to a periodic boundary but outside the 
computational domain Is obtained from a cell which also lies adjacent to a 
periodic boundary but Is Inside the computational domain. 

Since the flux Is zero on solid surfaces only the pressure need be known. 
This can be extrapolated from the Interior or determined from an adaptation to 
the present system of equations of a normal pressure gradient condition devel- 
oped by Rlzzl (ref. 6). The present work uses the adapted Rlzzl condition on 
the hub and extrapolation for the blade surfaces. 


Mesh Generation 

To solve the average-passage equation system through a multiblade row 
machine, a mesh Is needed for each blade row which contains the axial and 
radial coordinates of all blade rows. Thus, for a two-stage machine, four 
grids would be generated and each assigned the thickness and period of one of 
the four blade rows. However, each mesh must also conform axl symmetrically to 
the coordinates of the other blade rows. To do this efficiently, the geometry 
Is separated Into blade and nonblade sections from Inlet to exit. An 
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axlsymmetrlc algebraic mesh Is generated using two-dimensional spline fits and 
the axial and radial coordinates are common to all the grids In the meridional 
plane. To complete the grids, the tangential mesh lines are generated using 
spline fits and taking Into account blade thickness and blade count. 


Solution Procedure 

A nested Iteration procedure using an Inner and outer loop was presented 
In reference 9 to solve the average-passage equation system through each blade 
row of a multiblade row machine. Within the Inner loop the three-dimensional 
"average" flow variables are evaluated for a given distribution of body forces, 
energy sources, and correlations. These terms are denoted as S In equation (5). 
The Inner loop uses the Runge-Kutta Integration procedure outlined earlier for 
solving the Euler-like equations. The outer loop updates the body forces, 
energy sources, and correlations based on the axlsymmetrlc average of the con- 
verged Inner loop solution. An outline for updating the above terms can be 
found In reference 9. Global or outer loop convergence Is obtained when the 
differences between the axlsymmetrlc average of the time-average flow 
variables on each blade row Is below a given tolerance. 


RESULTS AND DISCUSSION 

A General Electric counter-rotating unducted fan configuration was sim- 
ulated on a Cray 1-S. The geometry contains two, elght-bladed fans (fig. 2) 
designed to operate at a free-stream Mach number of 0.72 and advance ratio of 
2.80. The grid (fig. 3) contains 99 axial, 36 radial, and 16 circumferential 
points with 28 points lying forward of the front blade, 20 points axially on 
both blades, 15 points between blades, and 15 points aft of the rear blade. 

Both blades contain 22 equally spaced points In the radial direction with the 
remaining 14 spaced from the blade tip to the free-stream. No mesh clustering 
was used In the angular direction. A sting whose diameter was equal to the 
sting at the hub exit was affixed to the front of the nacelle. 

Figures 4 and 5 show contours of constant relative Mach number on the 
pressure and suction side of the forward and aft blade rows. The shock loca- 
tion on the suction side of both blades (figs. 4(a) and 5(a)) run along the 
trailing edge of the blade terminating at about 30 percent of span on the first 
propeller. One would expect to observe this shock structure given the high 
Inlet Mach number and blade geometry. There also appears to be a shock In the 
vicinity of the nacelle/suctlon surface Interface for the forward propeller. 
Figure 5(b) shows the presence of a supersonic bubble at the junction of the 
nacelle surface on the pressure side of the aft blade. This bubble seems to 
extend across the passage to the suction surface Indicating that the flow In 
the hub region Is choked. 

The forward blade was also run In Isolation using the same spinner-nacelle 
geometry. The Inlet Mach number and rotational speed were Identical to those 
for the counter-rotating configuration. The relative Mach number contours are 
Illustrated In figure 6 for the suction and pressure surfaces. Note that the 
Mach numbers In general are lower than those which appear on the forward pro- 
peller of the counter-rotating configuration. This Is attributed to the 
"Induction" effect of the second blade row. More fluid passes through the 
counter-rotating configuration than through the Isolated configuration. This 
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can be seen more readily In plots of relative Mach number versus axial dis- 
tance. Three radial stations are shown In figures 7 to 9. It Is seen that the 
relative Mach number forward and aft of the first blade row Is higher for the 
counter-rotating configuration than It Is for the Isolated configuration. The 
shock strength for both configurations appear to be the same, however, for the 
counter-rotating configuration the shock location Is nearer the blade trailing 
edge. 


There are a number of ways to determine convergence of the above simula- 
tions. Jameson typically computes the time derivative of the density ( 1 . e . , 
3/o/at), and the number of supersonic points. Figures 10 show plots of 3p/3t 
versus the number of cycles for the first and second blade row simulation. The 
solid line Indicates the maximum absolute value of the derivative and the 
dashed line the average value of the derivative. The average value Is deter- 
mined by evaluating the sum of the absolute value of the time derivative at 
each point In the field divided by the number of points. The peaks at every 
500 cycles are attributed to updating the body forces, energy sources, and 
correlations. It Is seen that the present solution strategy for solving the 
Inner loop equations converges as Indicated by the reduction In both the maxi- 
mum and average levels of ap/dt. Figure 11 show the number of supersonic 
points based on the absolute Mach number versus the number of calculation 
cycles. Again, the number of supersonic points for both simulations converges 
to a constant value further Indicating convergence. Finally, figure 12 meas- 
ures the l _2 norm of the difference of the axl symmetrical ly-averaged solu- 
tions calculated at the end of each outer loop. The l _2 norm shows a drop 
of two orders of magnitude. This reduction was Judged sufficient to consider 
the computations converged. 

As a final check on the solutions described above, they were compared to 
experimental measurement taken In the NASA Lewis 8 by 6 Tunnel (ref. 10). This 
comparison Is shown In figure 13. The solid and dashed line Is a plot of the 
axlsymmetrlcally-averaged static pressure along the nacelle obtained from the 
first and second blade row simulations, respectively. The diamonds represent 
the experimental data. There Is good agreement between the prediction and 
measurement aft of the maximum nacelle diameter. The discrepancy at the 
spinner-nacelle nose was expected since the physical domain did not conform to 
the true nacelle geometry In this region. The grid conformed to a sting having 
the same diameter as the sting attachment at the end of the model. 

Since the above computations required large memory and CPU time, some code 
enhancements are being pursued. These Include multitasking and minimization 
of In-core storage. The first attempt at multitasking required the resources 
of a Cray X-MP multiprocessor computer. This effort Involved assigning the 
flow field simulation associated with each blade row to a processor. This step 
alone decreased run time by a factor of two for a two blade row configuration. 
Also, due to the advanced architecture of the Cray X-MP compared to the Cray 
1-S, an additional speedup was obtained. The net result was a reduction In CPU 
time by a factor of three over a comparable simulation on the Cray 1-S. To 
minimize In-core storage, the three-dimensional solutions will be stored In 
two-dimensional planes on secondary storage. . To minimize I/O overhead, the use 
of a high speed mass storage device will be needed. By Implementing these two 
enhancements, the capablllfles-of the code can be extended to solve multistage 
problems. 


8 



CONCLUSIONS 


A numerical procedure based on a finite volume formulation was developed 
to solve the average-passage equation system for a multiblade row configura- 
tion. This procedure employed a four-stage Runge-Kutta Integration scheme to 
march the equations forward In time towards the time asymptotic limit. A 
computer code based on this procedure was successfully used to simulate the 
average-passage flow fields associated with a high speed counter-rotating 
propeller. The results of this simulation yielded Information which proved 
useful In evaluating aerodynamic design. 
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Figure 2. - General Electric unducted fan configuration. 




(a) Constant Z, cut for forward and aft propellers. 



(b) Constant n cut for forward propeller. 



(c) Constant n cut for aft propeller. 

Figure 3. - GE unducted fan grid (99x36x16). 





(a) 10% span. 



(b) 50% span. 



(e) 90% span. 


Figure 7. - Relative Mach number vs axial distance on blade sur- ' 
face. Forward propeller. ■ : 




Zje-Zle 
(c) 90% span. 


Figure 8. - Relative Mach number vs axial distance on blade 
surface. Aft propeller. 
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Figure 9. - Relative Mach number vs axial distance on blade 
surface. Forward propeller in isolation. 
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(a) Forward propeller solution. 



NUMBED OF CYCLES 


(b) Aft propeller solution. 

Figure 10- Convergence history for forward 
and aft propeller. 


NUHBE P OF CYCLES 


(b) Aft propeller solution. 

Figure 11. - Convergence based on number 
of supersonic points. 
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